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Using a numerical implementation of strong disorder renormalization group, we study the low- 
energy, long-distance properties of layers and bilayers of S = 1/2 Heisenberg antiferromagnets 
with different type of disorder: bond randomness, site and dimer dilution. Generally the systems 
exhibit an ordered and a disordered phase separated by a phase boundary on which the static 
critical exponents appear to be independent of bond randomness in the strong disorder regime, 
while the dynamical exponent is a continuous function of the bond disorder strength. The low- 
energy fixed points of the off-critical phases are affected by the actual form of the disorder, and the 
disorder induced dynamical exponent depends on the disorder strength. As the strength of bond 
disorder is increased, there is a set of crossovers in the properties of the low-energy singularities. 
For weak disorder quantum fluctuations play the dominant role. For intermediate disorder non- 
localized disorder fluctuations are relevant, which become localized for even stronger bond disorder. 
We also present some quantum Monte Carlo simulations results to support the strong disorder 
renormalization approach. 



I. INTRODUCTION 

The two-dimensional (2d) spin-1/2 Heisenberg antifer- 
romagnet has attracted abiding interest in recent years, 
mainly motivated by its relation to high-temperature 
superconductivity^ According to the Mermin-Wagner 
theorem,— the Neel antiferromagnetic ( AF) long-range or- 
der in 2d can exist only at zero temperature, but even 
then it can still be reduced by quantum fluctuations. It 
has been established that at T = the AF order survives 
for several lattices, such as for the square lattice. The or- 
dered ground state is accompanied by gapless low-energy 
excitations, which, according to spin-wave theory^i and 
the non-linear cr-model description^, behave as: 

^E^^L~'% z, = 2, (1) 

where L is the linear size of the system, Zq is the dy- 
namical exponent and the subscript q refers to quantum 
fluctuations. The AF order in the ground state can be 
suppressed by introducing frustration (e.g. with diago- 
nal couplings in the square lattice: Ji — J2 model), ^ by 
dimerizing the lattice^ or by coupling two square lat- 
tices to form a bilayer /1^1^ By increasing these disorder- 
ing effects, the AF order is reduced progressively and will 
disappear at an order-disorder quantum phase transition 
point. 

In real materials impurities and other types of 
quenched disorder are inevitably present or can be con- 
trolled by doping. Fluctuations due to quenched dis- 
order can further destabilize the AF order, resulting 
in disordered ground states and random quantum crit- 
ical points. Quasi- two-dimensional materials, such as 



La2Cu04 doped with Mg (or Zn) and K2CUF4 (or 
K2MnF4) doped with Mg can be approximately de- 
scribed by the 2d AF Heisenberg model with static non- 
magnetic impurities. In these systems a disorder induced 
quantum phase transition from Neel order to a disordered 
spin liquid phase was observedii^ 

Theoretical investigations on the disorder effects in 2d 
Heisenberg antiferromagnets have been mainly restricted 
to dilution effects. Quantum Monte Carlo (QMC) sim- 
ulations of the diluted square lattice model showed that 
the AF long-range order persists up to the classical per- 
colation point and the critical exponents are identical to 
those of classical percolation for all S".^^ In studies of the 
square lattice model with staggered dimers and dimer 
dilution, unusual critical properties were found, among 
others, at the classical (bond) percolation point there is 
a critical line with varying exponentsii^ In the 2d bilayer 
Heisenberg antiferromagnet the random dimer dilution 
can be introduced by randomly removing the inter-layer 
bonds. In recent QMC simulations r^*i^*i^ random quan- 
tum critical points with an universal dynamical expo- 
nent z « 1.3 were deduced by varying the ratio of the 
inter-layer and intra-layer couplings below the percola- 
tion threshold. 

In the presence of bond randomness, the low-energy 
properties of the above mentioned 2d random models can 
be studied by a strong disorder renormalization group 
(RG) approachjii which was originally introduced by 
Ma, Dasgupta and Hu^^ for the Id random AF Heisen- 
berg model. In a detailed analysis of this RG procedure 
Fisher— solved the RG equations for the Id model analyt- 
ically and showed that during renormalization the distri- 



bution of the couplings broadens without hmit, indicat- 
ing that the RG flow goes to an infinite-randomness fixed 
point.-^ Due to infinite randomness, approximations in 
the RG procedure are neghgible and the scahng behavior 
of the system - both in dynamical and static sense - is 
asymptotically exact. The ground state of the Id model, 
the so-called random singlet state,-*' consists of effective 
singlet pairs and the two spins in a given singlet pair can 
be arbitrarily far from each other. Renormalization of 
the Id model with enforced dimerization (with different 
probability distributions of the even and odd couplings) 
leads to a random dimer phase^^ which is a prototype of 
a quantum Griffiths phase. The singular properties of the 
Griffiths phase are controlled by a line of strong disorder 
fixed points; along this line, the disorder induced dynam- 
ical exponent z varies continuously with the strength of 
dimerization. The dynamical exponent, calculated by the 
RG method, is presumably asymptotically exact, how- 
ever the static behavior, such as the density profiles, are 
correct only up to the correlation length in the system. 

Variants of the strong disorder RG method have been 
applied for various Id and quasi-ld (spin ladders) random 
Heisenberg models. In Heisenberg models with mixed 
ferro- and antiferromagnetic couplingsSS, during renor- 
malization large spins are formed and the dynamical 
properties of these large-spin phases are different from 
the Griffith phases, for example the uniform magnetic 
susceptibility has a Curie-like low-temperature behavior. 
The strong disorder RG method for more complicated ge- 
ometries, such as in 2d, can only be implemented numeri- 
cally and the calculated dynamical exponent z is presum- 
ably approximative. However we expect that the qualita- 
tive form of the low-energy singularities is correctly pre- 
dicted by these investigations. In previous studiea^i 2d 
and 3d Heisenberg antiferromagnets with/without frus- 
tration in the presence of bond disorder were numerically 
studied for random coupling constants taken from the 
Gaussian or from the box-like distributions. In contrast 
to the Id case, no infinite disorder fixed point is observed. 
Non-frustrated models are shown to have a conventional 
Griffiths- like random fixed point, whereas the dynamical 
singularities of frustrated models are controlled by large- 
spin fixed points. 

In the present paper we extend previous investigations 
of 2d random Heisenberg models in different directions. 
First, we consider strong disorder represented by power- 
law distribution of the couplings and study systemati- 
cally the variation of the dynamical singularities with 
the strength of bond disorder. In particular, we are in- 
terested in the localization properties of the low-energy 
excitations. Second, we consider non-magnetic impuri- 
ties and study the combined effect of bond disorder and 
site dilution. Our third direction of study considers AF 
bilayers with bond disorder and randomly removed inter- 
layer dimers. Evidently, with vanishing inter-layer cou- 
pling this problem reduces to our second model. 

The paper is organized as follows. The models un- 
der investigation as well as their basic properties are 




FIG. 1: The diluted bilayer model. Solid circles represent 
spins and open circles indicate the removed dimers. Neigh- 
boring spins in each plane interact with the coupling J, and 
the interplane coupling is K. 



presented in Sec. E] The strong disorder RG method 
and the properties of the basic fixed points are shown 
in Sec. IIIII A description of the QMC stochastic series 
expansion method, which is used to support the strong 
disorder RG approach, is given in Sec. II VI Results on the 
critical properties as well as the Griffiths singularities of 
different disordered Heisenberg AF models are presented 
in Sec.lVland discussed in Sec. I VII 



II. MODELS AND PHASE DIAGRAMS 

We start with the definition of the most general model 
considered in this paper: the double-layer Heisenberg an- 
tiferromagnet with random dimer dilution (see Fig. ^ 
which is described by the Hamiltonian: 



n=l,2 (ij) 



z,n ' ^j.n 



A^ 



Ki^i S.j,i • Si,2 -(2) 



Here Si^„ is a spin-1/2 operator at site i of the 71-th square 
lattice layer. The antiferromagnetic planar (inter-layer) 
coupling constants J^j {Ki) are independently and iden- 
tically distributed random variables. The dimer dilution 
at site i is represented by the variable e^, which is e, ~ 
with probability p and e^ = 1 with probability 1 — p. 

To our knowledge, this model has so far only been 
studied without bond disorder, i.e. Ki = K Wi and 
Jij = J yi,j. The schematic phase diagram of this 
model at zero temperature in terms of the coupling ratio 
g = K/ J and dilution p is shown in the plane £) = 
in Fig. 121 The point at (g = 0, p = 0) corresponds to 
two uncoupled non-diluted square lattice AF Heisenberg 
model and exhibits AF long-range order in its ground 
statei2^ At p = a finite inter-plane coupling, g > 0, 
causes a tendency for neighboring spins in the adjacent 
layers to form singlets and the AF order is therefore re- 
duced. If the coupling ratio exceeds some critical value, 
g > Qc, the system will undergo a quantum phase transi- 
tion from an AF state to a disordered state. This T = 
order-disorder transition is expected to belong to the 
universality class of the 3d classical Heisenberg model 



H(g.O,0) 



disordered 



.B(g^P„,0) 




PD(0,Pp,D) 



FIG. 2: Schematic phase diagram of the dimer diluted bilayer 
Heisenberg antiferromagnet, as function of the coupling ratio 
g, the fraction of the removed inter-plane dimers p, and the 
strength of bond disorder D. The disordered phase and the 
AF ordered phase are separated by a critical surface, indicated 
by dashed lines, which is located ai p < pp and gdp, D) , where 
Pp is the site-percolation threshold. In the model without 
bond disorder, D = 0, there are two unstable fixed points, H 
and P, as well as a stable bilayer fixed point B. In the diluted 
single layer g = with bond disorder, the phase boundary 
is located at the percolation threshold with universal static 
and strong disorder dependent dynamical critical exponents, 
indicated by the line of fixed points PD. In the AF ordered 
phase the dynamical exponent z is determined by quantum 
fiuctuations for weak disorder (indicated by a grey region at 
<7 = 0), whereas 2 is _D dependent for strong disorder. 



according to the cr-model description by Chakravarty et 
alA Results of recent QMC simulations are in accordance 
with this conjecture and the critical ratio is calculated as: 
gc « 2.5220.2iS& 

Along the horizontal axis of Fig. 12 i.e. with g — 
(and -D = 0) we have two uncoupled site diluted Heisen- 
berg AF planes. Increasing dilution suppresses AF order 
progressively and according to QMC results the quan- 
tum phase transition takes place at the classical site- 
percolation threshold, ^^ pp = 0.407. Furthermore, the 
critical exponents arc those of the classical percolation 
transitioniH Now having both dilution, p > 0, and finite 
inter-layer coupling, g > 0, the phase boundary gc(p) 
is monotonously decreasing with increasing dilution^^ii^ 
However even at the percolation threshold there is a fi- 
nite critical coupling: gdPp) = Qp ~ 0.16, and this 
fixed-point, marked by B in Fig. |21 is found to control 
the phase transition between the ordered and disordered 
phases for p > and g > Giiii^ This fixed point is a con- 
ventional random fixed point with power-law dynamical 
scaling and universal exponents. "'^^ 

In this paper we extend the space of parameters by 



introducing bond disorder, such that the intra-layer and 
inter-layer couplings are independent and identically dis- 
tributed random variables taken from the distributions: 



<J) = 



p{K) 



J, 



-l/D 



K, 



D 

-l/D 



D 



J-'+'^'', for < J < J„,ax; (3) 
■i^-1+1/^, for < X < i^„,ax , 



respectively. Here D^ — var[lnJ] — var[lni^] mea- 
sures the strength of disorder (var[a;] stands for the 
variance of x) and the control parameter is defined as 
g = ii'max/'/max- Notc that an uniform distribution cor- 
responds to D = 1. In particular we are interested in the 
properties of the phase diagram, the singularities at the 
phase transitions as well as the form of disorder induced 
low-energy excitations in the different regions. 



III. THE STRONG DISORDER RG METHOD 
AND ITS FIXED POINTS 

The strong disorder RG methodpii is an important tool 
to study random quantum systems. Here we recapitulate 
the basic ingredients of the method used for the 2d ran- 
dom Heisenberg antiferromagnet. 

The RG proceeds by eliminating at each step a term 
in the Hamiltonian with the largest gap separating the 
ground state and the first excited state. This decima- 
tion process generates new effective couplings between 
the remaining sites which are calculated perturbatively. 
For a lattice with more complex structure than a single 
chain, such as the bilayer antiferromagnet, the renormal- 
ized Hamiltonian contains effective spins of arbitrary size 
with a complicated correlated network and has both an- 
tiferromagnetic and ferromagnetic (F) couplings. The 
RG procedure for this Hamiltonian thus consists of two 
types of decimation rules, one for singlet formation (for 
equal-size spins with an AF bond), and one for cluster 
formation (for all other cases). Further details of the RG 
procedure can be found in Ref. 22 24 28]. 

During rcnormalization the energy scale Q, which is set 
by the cutoff the energy gaps of the effective Hamiltonian, 
is gradually decreasing. In the vicinity of the low-energy 
fixed point fl* -^ 0, the low-energy tail of the distribution 
of the gaps for a large finite system of linear size L follows 
the relation: 



P(A,n,L) = L-P|§,if 



-§ 



(4) 

which defines the gap exponent lu. The energy-scale and 
the length-scale is related by 51 ~ L~^ with the disorder 
induced dynamical exponent z. Note that with the ini- 
tial power-law distribution of the couplings in Eq. © the 
initial gap exponent is given by ijjq = — 1 + 1/-D. At a 
conventional random fixed point, we have uj/luq — 0(1), 
while at a infinite-disorder fixed point the distribution 
of the effective gaps broadens without limit, indicating 



u) /ujQ —>■ oo. If the low-energy excitations are localized, 
than the gap distribution for a fixed A is proportional 
to the volume of the system: P{A,il,L) ^ L"^. From 
Eq.Q), we obtain in this case: 



(5) 



1 



here an exponent z' is defined. Note that at an infinite- 
disorder fixed point the dynamical exponent z is formally 
infinite. 

Another characteristic feature of the fixed point is 
the typical size of the effective cluster moment, Scs — 
l^jiSil, which is determined by the classical correla- 
tion of the spins in the ground state, and the positive 
(negative) sign corresponds to an F (AF) coupling. ScS 
is expected to scale as Sett ^ L'^'^ . There are two types 
of fixed points concerning the value of Q: In some models 
the decimated spin pairs are typically singlets or the size 
of the effective spins has a saturated value, which yields 
^ = in the low-energy limit; In some models, mainly 
with frustration, large effective spins are formed and if 
ferromagnetic and antiferromagnetic couplings are deci- 
mated uncorrelated one obtains^^ C, — 1/2. This state is 
called the large-spin phase. 

In the RG method static correlations can be mea- 
sured by considering the staggered ground-state correla- 
tion function C{r) between two spins at distance r. This 
is defined as 



C(r)^a, -(7y.,S,.S,) 



(6) 



where r]ij = (^—lyi+Vi+^j+Vj ^ and r is measured by 
1-norm distance (also known as Manhattan distance): 
r ~ Vij = \xi — Xj\ + \yi — yj\. This choice was 
made for computational convenience; in the limit r -^ 
cx) it yields the same asymptotic behavior of C{r) as 
the one calculated with the Euclidean distance r — 
\J{xi — XjY + (j/j ~ 2/i)^- III o^i' R-GI scheme, the corre- 
lations of spin pairs which form an effective spin at each 
RG stage, are calculated by 



(Si ■ Sj) = aifcQ;j7(Sfe -Sf ), 



(7) 



where aife^;) = (Si(j) • S^^,))/(S^^,) ) are the proportion- 
ality coefficients for each spin. We assume zero corre- 
lation between two spins that do not form an effective 
spin. After accumulating the correlations between all 
decimated spin pairs, we divide the correlation for a given 
distance r by Irl? , which corresponds to the number of 
pairs a 1-norm distance r apart. Here we note that the 
RG results for static correlations are expected to be valid 
only in the vicinity of a (static) critical point. Thus the 
calculated correlation functions for the Id problem are 
asymptotically correct only in the vicinity of the phase 
boundary. 

Within the RG study, thermodynamics can be under- 
stood by stopping the RG procedure when the energy 
scale, i.e. the cutoff of energy gaps £7 in our case, reaches 



the thermal energy at a given temperature Tm^t^ At this 
scale, almost all decimated spins are effectively frozen, 
while almost all remaining spins involve couplings which 
are much less than T and hence can be regarded as free. 
The magnetic susceptibility per spin is then mainly given 
by the Curie-contribution of those remaining spins and 
is given by 



xiT) 



1 



riT 



TL^^^ 



S^{S, + 1) , 



(8) 



where the summation runs over all clusters left at the 
given temperature T, and Si is the (effective) spin mo- 
ment. In the low-temperature limit the susceptibility 
generally behaves as a power-law: 



xiT) ^ T~ 



(9) 



If during renormalization there is no large-spin formation 
i.e. C = Oi then 6 = oj va the low T limit, whereas in the 
large-spin phase with C, = 1/2 there is a Curie-like de- 
pendence: 9 = 1. Singularities of the specific heat or the 
magnetization can be calculated similarly, see Ref . \V^ . 



IV. QUANTUM MONTE CARLO METHOD 

A. Description of the method 

Here we use the QMC stochastic series expansion 
(SSE) method within a directed loop framework intro- 
duced by Syljuasen and Sandvik in Ref. 1291 Starting 
with a general Heisenberg Hamiltonian with random ex- 
changes J(6), we can rewrite it as a sum over diagonal 
and off-diagonal operators 



iVt 



n 



E^(^) 



fc=l 



Hi, 



H. 



2,6 



(10) 



where h denotes a bond connecting a pair of interacting 
spins (i{b) , j (b)) , and Nt is the total number of bonds. 



i?i,b - C- S'f((,)S'J(fc) (11) 

is the diagonal part and the off-diagonal part is given by 



1, 



Him - T;[St(b)^ 



2^~i('b) j{b) 



^i(b)^](b)\^ 



(12) 



in the basis {\a)} = {|5f , 51, ..., S'£}}. The constant 
C which has been added to the diagonal part ensures 
that all non-vanishing matrix elements are positive. The 
SSE algorithm consists in Taylor expanding the parti- 
tion function Z = Tr{e~^^} up to a cutoff M which is 
adapted during the simulations in order to ensure that 
all the elements of order higher than AA in the expansion 
do not contribute. So 



Z 



EE 



I3'\M ~n)\ 
7W! 



M 



a 



n Ji^^)Ha^ 



(13) 



where Sm denotes a sequence of operator indices 

Sm = [ai,bi], [a2,b2],--[aM,bM] 



(14) 



with Gi — 1,2 corresponds to the type of operator (diago- 
nal or not) and hi = 1, 2, ...Ni, is the bond index. A Monte 
Carlo configuration is therefore defined by a state \a) and 
a sequence Sm- Of course, a given operator string does 
not contain A4 operators of type 1 or 2, but only n; so 
in order to keep constant the size of Sm: M. — n unit 
operators i^o.o = 1 have been inserted in the string, tak- 
ing into account all the possible ways of insertions. The 
starting point of a simulation is given by a random ini- 
tial state I a) and an operator string containing M unit 
operators [0, 0]i, ..., [0, 0]^. The first step is the diagonal 
update which consists in exchanging unit and diagonal 
operators at each position p [0, 0]p ^^ [1, bi\p in Sm with 
Metropolis acceptance probabilities 



J{b)Ni,(3{a{p) m^b a{p) 
-P[o,o]p^[i,b]p = min(l, — ), (15) 



-P[i,6]p^[o,o]p = min(l. 



M -n 



J{h)Nb(3{a{p) 



Hi, 



a(j>) 



).(16) 



During the "propagation" from p = 1 to p = A^, the 
"propagated" state 



l«b))~n^a.,b,la) 



(17) 



i=l 



is used and the number of non-unit operators n can varies 
at each index p. The following step is the off-diagonal up- 
date, also called the loop update, carried out at n fixed. 
Its purpose is to substitute [1, bi\p ^^ [2, bi]p in a non- 
local manner but in a cluster-type update. At the SU{2) 
AF point, the algorithm is deterministic because one can 
build all the loops in a sole way.^^ One MC step is com- 
posed by one diagonal and off-diagonal updates. Before 
starting the measurement of physical observables, one 
has to perform equilibration steps, notably necessary to 
adapt the cutoff M . 



B. Monte Carlo measurement issues 

The precise determination of physical observables us- 
ing QMC suffers obviously from statistical errors since 
the number of MC steps is finite. As we deal with dis- 
ordered spin systems, the sample to sample fluctuation 
is another source of errors. However one can use a rel- 
atively small number of MC steps for each sample (typ- 
ically ~ 100 at each temperature) since for the strong 
disorders considered here, the sample to sample variation 
produces larger error bars than statistical errors. Than 
we need to perform a disordered samples average over 
a significant number of realization: typically we use 10'^ 
samples. 



In order to study the low temperature properties, we 
use the /3-doubling strategy introduced by Sandviliii to 
accelerate the cooling of the system during a QMC simu- 
lation. Such a scheme is a very powerful tool because 
it allows to reach extremely low temperatures rather 
rapidly and reduces considerably equilibration times in 
the MC simulation. The procedure is quite simple to im- 
plement and its basic ingredient consists in carrying out 
simulations at successive inverse temperatures /3„ = 2", 
n = 0, 1, ..., Umax- Starting with a given sample at n = 
we perform a small number of equilibration steps N^q fol- 
lowed by Nm = 2Neq measurement steps. At the end of 
the measurement process, P is doubled (i.e. n — > n -I- 1) 
and in order to start with an "almost equilibrated" MC 
configuration, the starting sequence used is the previous 
Sm doubled, i.e., 

S2M = [ai,bi], ...[aM,bM][aM,bM], ■■■dO'iibi]. (18) 



V. NUMERICAL RESULTS 

In practice we started with a finite system of linear size 
L (up to L = 64) with periodic boundary conditions for 
each single layer, and performed the RG procedure un- 
til the last effective spins (or the last spin singlet). The 
static characteristics of the system, in particular in the 
vicinity of the phase boundaries, can be deduced from 
the average spin-spin correlation function. On the other 
hand, the form of the dynamical singularities can be ob- 
tained from the temperature dependence of the uniform 
susceptibility and from the distribution of the first en- 
ergy gaps corresponding to the energy scale of the last 
decimation step. From the histogram of the gaps we have 
extracted the gap exponent w and the dynamical expo- 
nent z, as discussed in Sec. IIIII Depending on the size 
of the system we have considered 1000 — 10000 disorder 
realizations. 

For the single layer we also compare the RG results 
with QMC simulations performed at finite temperature 
on 32 X 32 square lattices and averaged over 1000 random 
samples. 

In what follows, we present the phase diagram of 
the system and the properties of the different bond- 
randomness driven phase transitions. The dynamical 
properties of the ordered and disordered phases are dis- 
cussed afterwards. 



A. Phase diagram and critical properties 

Our main results are summarized in the schematic 
phase diagram of the system depicted in Fig. [21 It con- 
tains two phases: the ordered AF phase and the disor- 
dered paramagnetic phase. The phase transition between 
these two phases is controlled by several fixed points as 
shown in the phase diagram. The fixed points located at 
D = 0, denoted by H, B and P in Fig.|21 had already been 



TABLE I: Critical exponents at the fixed points of the bilayer 
Heisenberg antiferromagnet with random dimer dilution and 
bond disorder, see in Fig. |5| H: non-random bilayer (clas- 
sical 3d Heisenberg model); P: diluted single layer (classical 
2d percolation); B: dimer diluted bilayer; PD: diluted single 
layer with bond disorder. In the last rows critical exponents 
measured at two general points of the critical surface are pre- 
sented. 
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fixed point 


position (g,p,D) 


P/u 


z 


V 


jjSO 


(gc, 0,0) 


0.51 


1 


0.70 


p31 


(0,Pp,0) 


5/48 


91/48 


4/3 


gis 


{9p,Pp,0) 


0.56 


1.31 


1.16 


PD 


iO,Pp,D >0) 


0.50 


~3.2D 






(f.2, 0.33, 0.7) 


0.56 


1.36 






(7.5- f0"'',0.33,3) 


0.80 


5.13 





carefully studied by QMC simulationsiSiiiiiii^ The mea- 
sured critical exponents at these fixed points are shown 
in Table ^ along with the results for Z? > obtained 
from our study. 

We first consider the fixed points (PD) at the perco- 
lation threshold p = pp ioi g = 0. Fig. 13 shows the 
average spin-spin correlation function Cav(?') at g = 
for different dilution p and for strong bond randomness, 
D = 3{D — 10). From p < pp to p > pp the decay 
of Ca.v{r) in the log-log plot changes from a upward to 
a downward curvature, and at the percolation thresh- 
old a power-law decay is found, which implies that the 
position of the phase transition in the site-diluted sin- 
gle layer Heisenberg antiferromagnet is robust against 
strong bond disorder. In comparison to the case with- 
out bond disorder, the decay of the critical average cor- 
relation function is, however, faster. From the decay of 
Cav{r = L/2) ~ ^-2/3/1/ -f-Qj. djffgj-gjii; system sizes L, the 
decay exponent is found independent of the strength of 
the disorder for D > 3, and estimated as 2(i/v = 1.01(7), 
while 213/v = 0.21 for D = 0. This indicates that the 
percolating cluster is no longer ordered in the presence 
of strong bond randomness:^ Unlike the decay exponent 
of CavC?"), which is IJ-independent, the dynamical expo- 
nent z' obtained from the slope of the gap distribution 
is found to depend linearly on the strength of the disor- 
der in the large D region: z k, 3.2D, as shown in Fig. 0] 
For weak bond disorder D < 1, instead, we find that z' 
approaches to the value z — 91/48 for the D — case. 

Now we turn to the phase boundary for finite bilayer 
coupling g > 0. At a given p < pp and a fixed Z?, we 
calculated the average spin-spin correlations Cavii") for 
different values of the bilayer coupling g. As illustrated 
in Fig. El for p = 0.33, we find that the decay behavior 
of Cuvir) changes its characteristic from the one for the 
AF ordered phase to the one for the disordered phase as 
a critical value of gc is traversed. For weak bond dis- 
order D = 0.7, the critical coupling is located around 
gc — 1.2 and we note that the decay exponent of the crit- 
ical correlation, 2/3/:/ « 1.12, is approximately the same 
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FIG. 3: Log-log plot of the average spin-spin correlation 
function a.t g — measured for a L = 64 lattice with bond 
randomness D = 3 and D = 10 (inset) for different site di- 
lution p. The data are scaled to unity at r = 1. For p > pp 
the curves show downward curvature, indicating a faster de- 
cay than a power law characteristic to the disordered phase. 
At the percolation threshold p — Pp, the decay exponent is 
2/3/i/ ~ 1.01(7), which does not depend on D. For p < pp the 
curves bend upward, indicating a finite limiting value and a 
characteristic of the AF ordered phase. 




FIG. 4: Disorder dependence of the z' exponent at the per- 
colation threshold (at the line of fixed points PD in Table U 
in a log- log plot. The slope of the straight line is unity. 



as for D = 0. For strong bond disorder D = 3, the phase 
boundary shifts to a very small value of gc ~ 0.00075 
with the critical exponent 2/3/z/ « 1.6. The extreme 
small value of gc, which decreases even with D, makes the 
investigation on Z3-dependence of the critical exponents 
difficult. From our results for Cav{r) up to 13 = 5, the 
decay exponent P/v appears to be Z3-independent for a 
given p in the strong disorder regime, while it varies with 
the dilution p. To locate the critical bilayer coupling gc 
we also made use of the results for the dynamical expo- 
nent z', cf. Fig. El for p = 0,p = 0.125 and 0.33. As 
g is increased, the dynamical exponent is approximately 
independent of the value of g, but jumps to another g- 
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FIG. 5: The in-plane average spin-spin correlations of the 
double-layer AF model versus r in log-log plots for bond ran- 
domness D = 0.7 (left) and D = 3 (right) for different values 
of the bilayer coupling at p — 0.33 for L = 48. The data 
are scaled to unity at r = 1. We observe a crossover from a 
upward curvature through a power-law decay to a downward 
curvature. The order-disorder transition point gc « 1.2 shows 
an asymptotically linear dependence in the large r regime with 
a slope 2/3/1/ « 1.12(4) (indicated by a broken line) which is 
approximately the same as for g = 0. For D = 3 the transi- 
tion point shifts to a very small value of gc = 0.00075 and the 
critical exponent is estimated as 2/3/;/ ~ 1.60(4). 



independent value around the transition point. For weak 
bond disorder, such as D = 0.7 for p — 0.33, we find 
z' « 1.36, which is close to the value found for the case 
without bond disorder—. For strong disorder, in which 
case the RG approach is expected to be more appropri- 
ate, the dynamical exponent increases with D, which is 
a tendency already noticed for g = 0. 

To summarize our numerical findings indicate two dif- 
ferent regimes of phase transition. For weak bond disorder 
the static critical exponent (i/v as well as the dynami- 
cal exponent z' seems to coincide with the values for the 
case without bond disorder. For strong bond disorder the 
critical coupling gc is reduced to a very small value, the 
static exponent approaches a D independent, but dilu- 
tion dependent value, whereas the dynamical exponent at 
the transition point depends (linearly) on the strength of 
the bond randomness. The position of the order-disorder 
transition for a single layer (corresponding to g = 0) is lo- 
cated at the percolation threshold. Along the line of PD 
fixed points, the exponent (3/v deviates from the value 
for D = 0, but seems to be U-independent, while the 
dynamical exponent exhibits a linear dependence on D 
in the large D limit. 



B. Griffiths singularities in the ordered phase 

As discussed in the preceding subsection, the random 
dimer diluted bilayer antiferromagnet exhibits AF order, 
provided p < pp and the bilayer coupling is sufficiently 
small. The low-energy fixed points governing the Grif- 
fiths singularities in the ordered phase are of different 
types in the specific regions. These fixed points are in 
turn an effective singlet for p = and g — (single 
layer without site dilution), a large-spin fixed point for 
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FIG. 6: Variation of the gap exponent with the coupling 
ratio g for weaker (left) and stronger (right) bond disorder 
for different values of the dimer dilution. Note that in the 
ordered phase g < Qc a.s well as in the disordered phase g > 
Qc, z' is approximately independent of g. For weaker bond 
disorder there is a jump at the transition point g = gc and 
the dynamical exponent is close to the value z{gc) ~ 1.3 at 
the fixed point B for D — case, which is denoted by a broken 
line. For strong disorder (right) the transition point is located 
at a very small g so that it cannot identified in the figure. 



< p < Pp and g = (single layer with site dilution), 
and an effective singlet for < p < pc and < g < gc 
(bilayer with dimer dilution). In the following we study 
these different cases separately. 



1. Two-dimensional undoped antiferromagnet 

We start by discussing the results for the two- 
dimensional random Heisenberg model, which corre- 
sponds to g = and p = 0. A recent numerical studyS 
suggested that the AF order in this region vanishes only 
in the limit of infinite bond randomness. In our prelimi- 
nary study24 we showed that the low-energy fixed point 
of the model is conventional, however, the dependence on 
the strength of disorder was not investigated extensively. 
Here we calculate the gap exponent uj and the related 
exponent z' defined in Eq.(|SJ), as well as the dynamical 
exponent z, as a function of the disorder strength D. The 
gap exponent oj is obtained from the slope of the distri- 
bution of the log-gaps in the small gap limit, whereas 
the dynamical exponent is determined from the optimal 
scaling collapse of the curves according to Eq.Q as il- 
lustrated in Fig. For localized excitations the scaling 
curve is conjectured^i from extreme- value statistics to be 
described by the Frechet distribution^ 



P,{u) 



= ^„,''/--i 



exp(-u''/^) 



(19) 



with d = 2 and u = uqL^A, where uq is a non-universal 
constant. 

Both z and z' have an approximately linear in- 
dependence in the strong disorder region (D > 3) as 
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FIG. 7: A scaling plot of the log-energy gaps for the 2d 
antiferromagnet with strong bond randomness D = 8 ob- 
tained from 10000 samples for each size. The gap exponent, 
uj « —0.64, follows from the slope at small energy gaps and 
the dynamical exponent, z ~ 5.5, is determined by the fit pa- 
rameter in Eq.(|3J. Note that the relation in Eq.Q is satisfied, 
implying that the low energy excitations are localized. The 
full line represents the Frechet distribution given in Eg. 1191 1. 



shown in Fig.|Sl while no significant disorder dependence 
{uj w 0.7) is found for weak disorder i^ The exponents z 
and z' are found identical only for quite strong disorder 
D > 7. This indicates that the low-energy excitations 
are localized only in the strong disorder regime. 

We note that the vanishing energy gaps calculated by 
the RG approach are solely induced by disorder. How- 
ever, quantum fluctuations also induce vanishing gaps 
which are characterized by a dynamical exponent, Zq = 2, 
see in Eq.fQ). The true dynamical exponent is then given 
by ztTue = max{zq,z}, so that ztruc = Zq = 2 for weak 
randomness D < 3. 

We have also calculated the uniform magnetic suscep- 
tibility as a function of the temperature, which is shown 
in Fig. O for different disorder strength. Both RG and 
QMC results are shown and they display an excellent 
agreement. For strong bond randomness, the low-T sus- 
ceptibility exhibits power-law temperature dependence 
as given in Eq.® and the exponent is disorder depen- 
dent (see table EJ . The same behavior of the magnetic 
susceptibility has been found for the antiferromagnetic 
spin-1/2 ladders?^ Note however that the QMC results, 
shown on the right panel of Fig. O display a slow satura- 
tion of X when T — > (at least for D < 5) which is not a 
finite size effect^ but a signature of a tendency towards 
Neel ordering at T = Oi2^ 

Finally, we note that effective spins with size larger 
than 1/2 are formed during RG procedure due to the gen- 
eration of F couplings. In the low-energy limit, the over- 
all strength of the F couplings, however, becomes much 
weaker than that of the AF couplings, which leads to 
the disappearance of large effective spins and the singlet 
ground state. This agrees with the Marshall's theorem-- 
which states that the ground state of a bipartite AF 
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FIG. 8: Variation of the disorder induced dynamical ex- 
ponent z and the exponent z' with the bond randomness 
strength D at g — and p = in a log-log plot. Note that the 
dependence for D > 3 is approximately linear and the values 
of z and z' fit well for D > 7, indicating that the low-energy 
excitations are localized. 



TABLE II: Exponent 9 of the divergence of the uniform sus- 
ceptibility for various disorder strengths D for p = g = 0. 
Comparison between RG and QMC estimates. 



D 


Org 


6qmc 


3 


0.36 


0.37 


5 


0.60 


0.61 


10 


0.77 


0.81 



Hamiltonian with equal size sublattices is a total spin 
singlet. 



2. Two-dimensional antiferromagnet with site dilution 

The low-energy behavior of the site-diluted Heisenberg 
antiferromagnet is controlled by a large-spin fixed point, 
which is different from the undoped case where the last 
decimated pair of spins is an effective singlet. The sit- 
uation is similar to that of antiferromagnetic spin-1/2 
ladders with random site dilution. In this case Sigrist et 
al'^^ argued that if two vacancies are in the same sublat- 
tice, the ground state is no longer a singlet thus there 
are effective spins of size larger than 1/2. This has been 
verified by numerical strong-disorder RG calculations!^ 
In the 2d site-diluted case we also observed in our numer- 
ical RG calculation that the energy gap associated with 
an effective F coupling may become the largest gap to 
be decimated at some stage of the RG, especially in the 
low-energy regime. This will then lead to the formation 
of large effective spins as described in Sec. Illll 

We calculated the average size of the effective spin at 
the last decimation step (Scs) for various dilution con- 
centrations {p = 0.125, 0.33 and at Pp) and system sizes 
L. In the ordered phase, below the percolation threshold. 
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FIG. 9: Disorder average uniform susceptibility xiT) as a function of temperature T for various disorder strength D at 
g = p — 0. Left panel: RG results. From the low-temperature regime (T < 10~^) the exponent 6 is estimated as: 9 — 0.36 for 
D = 3, = 0.60 for D — 5, 9 = 0.71 for D = 7 and 9 — 0.77 for D — 10. For all cases studied, the temperature dependence 
deviates from Curie-like 1/T behavior indicated by the dashed line. Right panel: QMC results obtained on systems of 32 x 32 
spins. The exponent 9 is estimated in a range of T € [T* , 1] as: 9 = 0.37 for D = 3 (T* ~ 0.02), 9 = 0.45 for D = 3.5 
(T* ~ 0.01), 9 = 0.52 for D = 4 (T* ~ 0.01), 9 = 0.61 for D = 5 (T* ~ 0.002) and 9 = 0.81 for D^W (T* ~ 0.0001). 
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FIG. 10: Variation of the disorder averaged spin size (Scb) 
with the linear system size L in a log-log plot for different 
site dilution p at D = 3. for the single layer antiferromagnet 
g = 0. For p < pp the spin size follows (Scft) ~ L, indicated 
by the broken lines, whereas at p = pp the asymptotic power 
for large system sizes agrees with 91/96. 




FIG. 11: Temperature dependence of the uniform suscep- 
tibility per size for a diluted single layer, <; = 0, in log-log 
plots for various dilution concentrations p and for different 
bond random strength D = 3 (left) and D = 10 (right). The 
Curie- like 1/T behavior is indicated by straight lines. Both 
RG (upper panels) and QMC (lower panels) are shown 



P ^ Ppj the average spin size is found to increase linearly 
with the system size: 



('5'cff) ^ L, 



(20) 



which is demonstrated in Fig. ^] This result agrees with 
the scenario for the large-spin phase, as discussed below 
Eq.©. At the percolation threshold the same argument 
leads to (5'off) ^ L'^f/^, with df = 91/48 being the fractal 
dimension of the percolation cluster^^. 

A hallmark of the large spin phase is the universal tem- 
perature dependence of some thermodynamic quantities, 
in particular the disorder averaged uniform susceptibil- 



ity given in Eq.@ shows a Curie- like behavior at low 
temperatures. This is checked in Fig. 1111 in which the 
susceptibility obtained from both RG and QMC is plot- 
ted for different strength of bond randomness and dilu- 
tion concentrations. For not too strong bond disorder the 
agreement with the Curie-law is good, while for strong 
bond disorder this agreement is observed only for very 
low temperatures. Note that in the undoped regime the 
susceptibility exponent is a continuous function of the 
disorder, see in FigEI 

From the distributions of the low-lying energy gaps. 
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FIG. 12: The gap exponent uj in the diluted single layer 
antiferromagnet {g = 0) for different dilution and bond disor- 
der. Note that in the disordered phase, above the percolation 
threshold p > Pp, the gap exponent is practically independent 
of the dilution. 



we obtained the dynamical exponent z and the gap ex- 
ponent to. Unhke for the undoped model, the exponent 
z and z' = 2/(1 -I- u) in general do not agree with each 
other for < p < Pp, even in the regime of strong bond 
disorder. This indicates that fow-energy excitations are 
not localized due to the formation of large spins. Fig. 1121 
presents the exponent z' as a function of p for D = 3 
and D = 10. For a given bond disorder, the exponent 
varies continuously with p in the ordered phase (p < Pp), 
while going approximately to a constant in the disordered 
phase {p > Pp). We remind that to obtain the true dy- 
namical exponent ztruo one should also consider the effect 
of quantum fluctuations and thus ztruc = max{zq,z} in 
the ordered phase. 



3. The double-layer Hetsenberg antiferromagnet 

In the presence of random bilayer couplings g > 0, the 
low energy properties of the ordered phase are controlled 
by an effective singlet, both for p ~ and < p < Pc, 
which in turn is the same as for the single layer undoped 
model, see in Sec lVBTl Indeed, we observed similar 
low-energy properties. The dynamical exponent z, and 
the exponent z' are disorder dependent, but vary only 
weakly with the bilayer coupling g, see FigEl z and 
z' are identical only for strong enough disorder, when 
the low-energy excitations are expected to be localized. 
The average uniform susceptibility has a disorder depen- 
dent low-temperature behavior and the exponent 0, cor- 
responds to the gap exponent w. 



C. Griffiths singularities in the disordered phase 

The disordered phase of the system is divided into two 
parts with different low-energy properties: 
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FIG. 13: A finite size scaling plot of the distribution of the 
logarithm of the energy gap for the double-layer antiferromag- 
net with a bilayer coupling g = 1, bond randomness D — 8 
and dimer dilution concentration p — 0.125. The dynamical 
exponent z and the slop ( — 1 — uj) of small energy gaps agree 
well with the relation z = 2/(1 + 0;), implying localized energy 
gaps. 



• Above the percolation threshold p > Pp and g — 
the spins form only finite connected clusters. As a 
consequence the average effective spin has a finite value, 
as shown in Fig. ^|for p = 0.42. Due to the unpaired 
spins in the isolated connected spin clusters the average 
uniform susceptibility is Curie-like (see Fig. ^] for p = 
0.5). The dynamical exponent z and the gap exponent uj 
depend approximately linearly on the bond disorder D, 
they exhibit however no significant dependence on p, as 
shown in Fig. [T^ 

• Above the critical bilayer coupling g > gdp, D), the 
ground state is an effective singlet and in accordance with 
this the low-temperature uniform susceptibility is char- 
acterized by a non-universal exponent 6. For p < pp, 
there is a infinite cluster and the low-energy physics is 
governed by rare finite regions which are locally ordered. 
The low-energy excitations connected to these regions are 
thus expected to be localized, provided the bond disor- 
der is sufficiently strong. This is illustrated in Figll3l 
in which the scaling collapse of the energy gap distribu- 
tion is obtained for z = z' in accordance with Eq. I^J . For 
p > Pp, the connected spin clusters are finite and isolated. 
Therefore the low-energy excitations are also localized. 



VI. SUMMARY AND DISCUSSION 

In this paper we have studied the effect of strong 
bond disorder on the low-energy, long-distance proper- 
ties of Heisenberg antiferromagnetic layers and bilayers 
with site and dimer dilution. In particular we are inter- 
ested in the structure of the phase diagram, the form of 
the critical singularities as well as the properties of the 
Griffiths singularities. 
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In a single layer the order-disorder transition point is 
found to be at the percolation threshold p = Pp, thus for 
p < Pp the AF order survives for any finite value of bond 
disorder strength D. In contrast to this, the AF order 
at the percolation cluster, which is present for D = 0, is 
destroyed by bond disorder. This information is deduced 
from the decay of the average spin-spin correlation func- 
tion which has the same power-law form with a strong- Z? 
independent exponent "ifi/v, but much smaller than the 
known exponent 2l3/v = 10/48 at D = 0. The dynami- 
cal exponent z of the diluted single layer is found to be a 
continuously increasing function of the disorder D. Here 
we note that in the limit of infinite D the fixed point be- 
comes an infinite disorder fixed point with z — > oo so that 
the RG method is expected to be asymptotically exact 
with increasing D, as well supported by comparing with 
QMC results. 

In the dimcr diluted bilayer with g > 0, weak disorder 
is found not to modify the static critical exponent P/i' 
as well as the dynamical exponent z, which are - within 
the error bars - the same as one measures at the fixed 
point B. On the other hand for strong bond disorder 
the critical bilayer coupling is reduced to a very small 
value and both the static and the dynamical exponents 
are different than for weak disorder. While the static 



exponent approaches a D independent limiting value the 
dynamical exponent shows a linear D dependence. 

Considering the Griffiths singularities the low-energy 
fixed point of the RG is found to depend on the specific 
form of the disorder. For example, the non-diluted single 
layer {g = p = 0) transforms into an effective singlet, the 
diluted single layer {g = 0, < p < Pp) into a large spin, 
whereas the dimer diluted bilayer also into an effective 
singlet. In each cases the disorder induced dynamical 
exponent is found D dependent for sufficiently large D. 
For smaller values of D the true dynamical exponent is 
determined by quantum fluctuations, so that in this re- 
gion disorder can influence only the corrections to scaling. 
The low energy excitations are found to be non-localized 
for weak bond disorder as well as in the large-spin phase, 
and become localized only for substantially large disor- 
der. 
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